%This file generates a case of a 20x20 Cov-matrix with up to 10 sub blocks.
%sampled sequences with different levels of noice are generated and used 
%as estimated subblocks. The performance is measured for different number 
%of blocks and different noise levels. Sample length is fixed to 5000.

clear, close all
NumOfStocks = 20;
NumOfSamples=300;
MaxNoiseLevel=30;

%[CovBlocks, Combinations, SampleSizes, Sigma, StockReturns]=GetCovBlocks(0);
figure(1)
load data_10_blocks.mat

%Create a Sampled sequence of numbers. Maximum sample number = 10000.
Samples=Sigma^0.5*randn(NumOfStocks,NumOfSamples);
SampleNoise=randn(NumOfStocks,NumOfSamples);
error=[];
errorNoise=[];

for iterBlocks=2:5

    NumOfBlocks=2*iterBlocks;
    %Define number of samples per block
    SampleLengths=NumOfSamples*ones(1,NumOfBlocks);
    %Define the covariance blocks for this iteration instance
    CovBlocksIter=CovBlocks(1:NumOfBlocks);
    CombinationsIter=Combinations(1:NumOfBlocks);
    
    %Define weights according to sample length, make them sum 1.
    w=zeros(1,NumOfBlocks);
    for i=1:NumOfBlocks
        w(1,i)=SampleLengths(1,i)/sum(SampleLengths);
    end
    
    %Create the sampled Covariance blocks without noise
    CovBlocksSampled = CovBlocksIter;
    for k=1:NumOfBlocks
        CovBlocksSampled{k}=cov(Samples(CombinationsIter{1,k},1:SampleLengths(1,k))');
    end


    %Generate the optimal estimate when no noise is added
     cvx_begin
        cvx_quiet(true); 
        variable Sigma_hat(NumOfStocks, NumOfStocks) symmetric;

        %Define objective function
        f = 0;
        s_max = 0;
        for q=1:NumOfBlocks
            f = f + w(1,q)*norm(Sigma_hat(CombinationsIter{q}, CombinationsIter{q})-CovBlocksSampled{q},'fro');           
        end
        minimize (f)

        subject to

        Sigma_hat == semidefinite(NumOfStocks)
    cvx_end
    error(iterBlocks-1)=norm(Sigma-Sigma_hat,'fro');
        
    index=1;

    %Loop over different noise levels
    for iterNoise=0:10
        
        %Disturb the samples with noise
        SamplesIter=Samples+iterNoise*0.1*eye(NumOfStocks,NumOfStocks)*SampleNoise;
                
        %Create the noisy Covariance blocks
        CovBlocksSampledNoisy = CovBlocksIter;
        for k=1:NumOfBlocks
            CovBlocksSampledNoisy{k}=cov(SamplesIter(CombinationsIter{1,k},1:SampleLengths(1,k))');
        end

            cvx_begin
                cvx_quiet(true); 
                variable Sigma_hatWithNoise(NumOfStocks, NumOfStocks) symmetric;

                %Define objective function
                f = 0;
                s_max = 0;
                for q=1:NumOfBlocks
                    f = f + w(1,q)*norm(Sigma_hatWithNoise(CombinationsIter{q}, CombinationsIter{q})-CovBlocksSampledNoisy{q},'fro');           
                end
                minimize (f)

                subject to

                Sigma_hatWithNoise == semidefinite(NumOfStocks)
             cvx_end
             


        %Now test the result using frobenius norm

        errorNoise(index,iterBlocks-1)=norm(Sigma-Sigma_hatWithNoise,'fro');
        index=index+1;
        iterNoise
    end

    AbsoluteError(:,iterBlocks-1)=errorNoise(:,iterBlocks-1) - error(iterBlocks-1);
end

figure(1)
    hold on
    
    plot([0:10],errorNoise(:,1),':k')
    plot([0:10],errorNoise(:,2),'-*k')
    plot([0:10],errorNoise(:,3),'--k')
    plot([0:10],errorNoise(:,4),'-.k')

    legend('4 Blocks', '6 Blocks', '8 Blocks', '10 Blocks')
    title('Estimation error')
    xlabel('Level of noise on submatrices')
    ylabel('norm(S-Shat)')

figure(2)
    hold on
    
    plot([0:10],AbsoluteError(:,1),':k')
    plot([0:10],AbsoluteError(:,2),'--ok')
    plot([0:10],AbsoluteError(:,3),'--k')
    plot([0:10],AbsoluteError(:,4),'-.k')

    legend('4 Blocks', '6 Blocks', '8 Blocks', '10 Blocks')
    title('Increase in estimation error due to noise')
    xlabel('Level of noise on submatrices')
    ylabel('increase in error')
